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Abstract: 



A numerical solution is presented for the non-linear evolution equation that governs the 
^ ' dynamics of high parton density QCD. It is shown that the solution falls off as e^**/^ at large 

■ values of the impact parameter b. The power-like tail of the amplitude appears in impact 

parameter distributions only after the inclusion of dipoles of size larger than the target, a 
configuration for which the non-linear equation is not valid. The value, energy and impact 
parameter of the saturation scale Qs{y = ln(l/x), b)) are calculated both for fixed and running 
QCD coupling cases. It is shown that the solution exhibits geometrical scaling behaviour. The 
radius of interaction increases as the rapidity in accordance with the Froissart theorem. The 
solution we obtain differs from previous attempts, where an anzatz for b behaviour was made. 
The solutions for running and fixed as differ. For running we obtain a larger radius of 
interaction ( approximately twice as large), a steeper rapidity dependence, and a larger value 
of the saturation scale. 
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1 Introduction 



The main goal of this paper is to find a numerical solution (including its impact parameter 
dependence) of the non-linear evolution equation P], that governs the dynamics of the dipole 
scattering amplitude in the saturation region O El El • 

Although, there are several schemes for finding a numerical solution El 13 E] , we still lack 
a reliable method that will allow us to study the properties of the solution of the non-linear 
evolution equation. For example, there is the difficulty of understanding the large impact 
parameter [b = \b\) behavior of the solution. Our notation is such that the symbol "— 
represents a two dimensional vector, in the transverse plane. 

The large b behavior is strongly affected by the non-perturbative contributions 9], since in 
pQCD the amplitude falls off as a power of b. A decrease of this type leads to a power-like 
growth of the interaction radius, which violates the Froissart bound ^Uj. This difficulty can 
be overcome by introducing non-perturbative corrections for large b ( > l/m.,^), where m-n- 
is the mass of the lightest hadron (vr-meson). Two suggestions of how to include such non- 
perturbative corrections have been made : (i) to change the kernel of linear and non-linear 
equations [Hj; and (ii) to accumulate all non-perturbative corrections in the initial conditions 
of the equation [TTl|T2llI3j- The first attempt to solve the full non-linear equation including 
the impact parameter behavior, was made in Ref. The solution of Ref. ^3], shows that to 
obtain the correct behavior at large values of b, one has to alter the kernel of the equation. 

We examine this solution and show that the entire power-like decrease obtained in Ref. ^11 , 
originates in the kinematic region where the non-linear equation is not trustworthy. 

In our search for the solution we are guided by: (i) the known analytic expressions for the 
simplified cases Uni CSl , and (ii) practical experience gained in previous attempts: 

1. The non-linear evolution equation can be viewed as the sum of the semi-enhanced, 
'fan' diagrams (see Fig.^ for such a diagram). However, we should restrict ourselves 
to summing such diagrams only for the scattering amplitude of a dipole of transverse 
size (1?^), which is much smaller than the size of the target (21 El d'^ -C -R, where R is 
the target transverse size)^. If the target is smaller than the projectile, we need to sum 
diagrams of the type shown in Fig.^. Although, the kernel of the non-linear equation 
does not depend explicitly on the size of the target, this quantity appears in the solution 
by way of the region of validity of the equation. 

2. In the saturation region (r < 2/Qs, where Qs is the saturation scale) the solution dis- 
plays geometrical scaling [121 [1211211 Ej , namely, the scattering amplitude ^ is a function 
of only one variable 



We therefore expect that the entire b dependence can be absorbed into the saturation 
scale Qs{x; b). 

^For more detailed discussion of this key property of the non-linear equation see Refs. |19[ I20| . 
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3. We examine the behaviour of the scattering amphtude yf, b; xj, both for the case of 
the running QCD couphng constant, and for fixed a^. 

4. In the saturation region N (f,b;x^ < 1, consequently the non-hnear equation has a 
solution that satisfies s-channel unitarity. On the other hand, the behavior of the total 
cross section: 

a<^'P°i^(f,x) =21 £bN{f,bt;x) (2) 

depends critically on the large b dependence of the amplitude, and has to be studied 
separately. 

5. The semiclassical solution ^H] to the non-linear evolution equation, indicates that the 
behaviour at large b stems from the initial conditions. 

Our strategy is to assume that the size of the interacting dipole is smaller than that of the 
target. In the next section we describe the method of solution which we develop based on this 
strategy. In Section 3 we present our results. In particular, we check whether the geometrical 
scaling behavior in the form of Eq. (jT)) occurs in our solution. In the Conclusions, we discuss 
the main results and the possibility of using them for phenomenology. 

2 Method of the Solution 

The nonlinear evolution equation pP characterizes the low x behavior of the parton densities, 
while taking into account hdQCD effects, so that unitarity constraints are inherently obeyed. 

We denote by N the imaginary part of the interaction amplitude of a target with a parent 
dipole, of transverse size r and two dipoles, of transverse sizes ri and r2 {r = ri + r2), produced 
by r . 

The probability for the decay r ^ ri,r2 is given by the square of the wave function of the 
parent dipole, which, in a simplified form, can be written as r'^/r\r\. Each of the produced 
dipoles can interact with the target independently, with respective amplitudes iV(?/, b— \ri) 
and N{y,ri;b — ^^2), where y denotes the rapidity variable {y = — Inx) and b the impact 
parameter in the transverse plane. Adding these contributions, clearly overestimates the dipole- 
nucleon interaction, since there is the probability that during the interaction, one dipole can be 
in the shadow of the other. This correction is given by an additional negative quadratic term 
-N{y,r2;b- ^fi)N{y,fi;b - irg). 

In the procedure leading to the construction of the non-linear evolution equation, only fan 
diagrams of the type shown in Fig.^ were taken into account, while fan diagrams having 
multiple emissions from the parent dipole (Fig.H^) were omitted PP. In other words, the 
equation does not include contributions for |ri|, |r2| > R, where R is the transverse size of 
the target. The solutions, that have been presented previously in Refs. 001111111, included 
contributions from both large and small size dipoles'''. 

■fWe would like to mention that the point of view presented in Ref. [22] is very similar to ours. 
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Figure 1: Multiple emissions of small x gluons: (a) from the target, and (b) from the dipole. 



Our approach is initially to exclude contributions from large dipoles, and then to estimate 
the accuracy of this proceedure by employing a second iteration. The price we pay for this 
approach is, that our calculations are only valid for short distances, compared to the size of the 
target. Nevertheless, as we shall demonstrate later, contributions coming from large dipoles, 
which we consider as a correction, are relatively small. 

Hence, we can write the evolution equation in a differential form as: 

Ny{y,r;b) = — = j d r2^J{R - 9{R - Ir^]) 

2N{y,n-b - ^r^) - N{y,ri;b - lr2) N{y,r2]b - ^n) - N{y,r-b)] . (3) 

The r2 integration in Q is regulated by introducing an ultraviolet cutoff, which does not affect 
the physical quantities. 

Strictly speaking, the evolution equation is defined for fixed strong coupling, as- Since for 
running as it is difficult to define the running scale in terms of the variables of integration, fi 
and r2, of the r.h.s of Eq. Q. In the semiclassical approach, however, one sets as = as(4/r^), 
i.e. one uses the size of the parent dipole, as the running scale of the coupling. In the following, 
we shall attempt to compare the results obtained with fixed coupling to those obtained with 
running coupling in the semiclassical approach. Another approach would be to calculate the 
solution using as = dsH^rir 2/(^1 + ^2))^), but such a calculation is beyond the scope of the 
present work. 

The strong coupling constant ag is chosen to run according to the prescription 

An 



Po ln(4/(A^^^r2)) 

where Aqcd = 0.2 GeV and Po = (33 - 2nf)/3. For the case of fixed coupling constant we take 
as = 0.2. 

To obtain a solution to (jH]) at arbitrary b, one needs to specify the initial conditions of at 
y = yo, the rapidity from where we commence the evolution. The initial condition is a function 



3 



of both r and h aX y = yo, while our goal is to obtain a solution, which is a function of r and b 
for all y > yo- We are aware that during the evolution process, as y increases, the role played 
by the initial input decreases, and the evolution becomes the dominant factor influencing the 
dynamics per se. 

The experience we have in solving Eq. (jS)) (both with fixed and running as) stems mostly 
from the solution of the equation, where shifts in b on the r.h.s. of the equation, were neglected. 
For such a simplified equation, the impact parameter is treated as a fixed parameter in the 
procedure. For a solution of the type found in 0, once a trial solution was found at sufficiently 
large y, the curve approximating the solution was rescaled at fixed b assuming geometrical 
scaling [see Eq. ((T))], so as to reproduce the initial conditions at y = yo- The reconstructed initial 
conditions were then used as input for the evolution, so that, the final solution almost completely 
decoupled from the initial input. It was demonstrated in ^ that the initial conditions obtained 
in this fashion, at very large distances, provide a smooth extrapolation of the Glauber formula 
to unity. 

In this paper we adopt the basic principle of this procedure, and use the rescaled solution 
of jnj to construct the initial conditions. At first sight, it appears sufficient to follow jHj by 
adopting the fixed b solution as input, and assuming the following ansatz for the 6-dependence: 

N{yo, r- 6) = 1 - e-''^^"'^^ ^('^^Z^^") , (4) 

where S{b) is the profile function {e.g., a Gaussian or a dipole-like profile) and 

K{yo, r) = - In (1 - iV(yo, r; 6 = 0)). (5) 

However, in the Born approximation, at a finite mass of the gluon field, the scattering amplitude 
for large size parent dipoles, decreases already at 6 = 0. To facilitate such behavior, we modify 
(jH) as follows: 

iV(yo, r; b) = l- e-^^y^'"^ s(Vy^)/s(o) (g) 

Figs. EK-b show a comparison for the initial conditions of Eqs. and (jH)). The replacement 
6^ — > b^ + r'^, suppresses the initial conditions at large r, thereby smoothing the 6 functions 
of®. 

The most troublesome contribution in Eq. comes from the region of ri and/or r2 of the 
order of 2b. Indeed, in this region at large values of b, the r.h.s. of the master equation can be 
rewritten in the form 

2 

= ^ / d'^b (2iV(y,6,A6) - N'iy,b,Ab) - N{y,r,b)) (7) 

One can see that AA^ oc 1/6'^ at large b. This contribution has been discussed in Ref. [HIE]- 
We claim that this contribution originates from the interaction of the dipoles of size much larger 
than the size of the target, and that it cannot be evaluated using Eq. Q. 

Apart from rapidity, the function A^ depends on the dipole's degrees of freedom, which 
are the transverse separation of the parent dipole, r, and the transverse impact parameter. 
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Figure 2: The solution of the non-linear evolution equation for different initial conditions: 
(a) using Eq. @ with a Gaussian profile function; (b) using Eq. Q with the same 
profile. 



b. We therefore, define the following variables, through which we express the dipole-nucleon 
amplitude: r = \f\, b = \b\, and f|| = b ■ f/{br). In the following, we shall concentrate on the 
integrated quantity N{y, r;b) = J df\\N{y, r;b;f\\). 

Once the initial conditions were established, we applied the following numerical procedure, 
which has four stages: 



1. First iteration: we denote by yi a particular rapidity, at which the solution and its 
derivative are known for all values of r , 6 and f\\. The solution at = yi + h{y) was 
constructed as a matrix, in which the matrix elements correspond to (fixed y) solutions 
at different r and b, integrated over f\\: 

N^'\r,y,+,;b) = N^'\r,yf,b) + h{y,) J dr\\ iV«(y, r; 6; f||)|^^^^ , (8) 

where A^^^^ is given by the R.H.S of ©, and h{yi) is a variable step size in the rapidity 
space. The first iteration included successive constructions of such matrices over a range 
of about 10 units of rapidity, starting from ^ 4.6. The value of h{yi) was chosen 
dynamically to ensure that the maximal error over the ith matrix is small (we set our 
accuracy condition to 10~^). We used the Euler two-step procedure^ to select h{y), we 
found that the optimal step size may vary from h ^ 0.01 for small rapidities to h ^ 0.2 



^In the Euler two-step procedure, one obtains a solution along the path yi — > yi + h and an additional solution 
along the path yi yi + h/2 yi + h. The difference per unit step between the solutions is proportional to 
h{y)d^ N / dy'^ , and the convergence criterion is that the maximal difference over the matrix is small. 
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for large rapidities. The CPU time for the entire iteration was about 20 hours for 30 x 30 
matrices. 



Our main claim is that the first iteration as described above, provides a reliable solution 
to the master equation. 

2. First Correction: in order to check the accuracy of the first iteration, and to estimate 
the importance of long distance contributions, we substituted A^^^-* back into the non-linear 
evolution equation, and recalculated Ny in the kinematical region which was omitted from 
Q, namely for \r2\ > R- In other words, we calculated the contribution given by 
Eq. ((Tj). We consider the result of this calculation an estimate of the correction to A^^^-*, 
and denote it by AA^^^-*. 

3. Second iteration: to estimate the influence of the neglected terms on the behavior of 
the solution in the saturation region, we calculated the second iteration. Namely, we used 
the calculated AN^^^ (see Eq. ([7j)) as a second iteration, by adding the correction to the 
r.h.s of dHl): 

Ar(2)(^,_,,,r;6) = Ar(2)(^„r;6) + %,)(AiV«(y,,r;6)+ / df|| A^f (y, r; 6; f||) | ^ ) . (9) 

4. Second correction: finally, we calculated ANjj^\ the correction to the second iteration, 
by performing the integration over long distances as before. 



This procedure gives us full control both on the accuracy of the solution, and on the region 
of its applicability. The corrections of the first and second iterations are shown in Fig. El for 
6 = and 6 = 2 GeV~^, as a function of the ratio r/R. For not too large rapidities {y < 10), 
the corrections are of order 10% , even when the transverse size of the dipole r, is comparable 
to the transverse size of the target, R. For larger values of y, the corrections are small only for 
short distances. The second iteration is important at large rapidities where |AA^^^)| < lAA^^-^-*!. 
On comparing AA^ for both iterations, we see that we can only trust our solution for r/R < 0.4 
at large values of rapidity {y > 10), while for small values of y the accuracy of our solution is 
satisfactory, even for values of r/R ^ 1.5. This observation is most encouraging as we would 
like to use the solution to describe experimental data, which are mostly concentrated at y < 10. 



3 Evolution Results 
3.1 r, y and b dependence: 

Fig. m shows the evolution process from y = 4.6 to ?/ = 14 at 6 = 0, 1, 3 GeV^^ as a function 
oi r/R, and Fig. El shows the energy dependence of the amplitude at 6 = 0, for different values 
of r/R. We have chosen a target size of R^ = 3.1GeV~^ so as to compare with [Hj, where 
the fixed b solution was successfully fitted to the HERA experimental F2 data, using R^ as 
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Figure 3: The corrections to the first and the second iterations for 6 = and 6 = 2 GeV ^ . 

a parameter. In the next section, we discuss in more detail the differences and similarities 
between the present work and the solutions of [Hf I14j. 

Referring to Fig.|3J we see that at the beginning of the evolution process, the amplitude 
increases mildly with the dipole size, and does not exhibit saturation in the region oi r < R. 
As y increases, the rise of becomes steeper. This rise is tamed when the impact parameter 
becomes large. 

In our approach, only the kinematical region r/R^l can be explored accurately , although, 
there are some other kinematical regions where the corrections AiV^^^ and ANjj^^ are relatively 
small. At long distances there remain considerable theoretical uncertainties. 

Fig. [7| shows the 6-dependence of our solution for r/R = 1. It was argued in ^3], that for 
small rapidities the evolution process starts to generate a power behavior of b~^'^, which becomes 
even milder at higher rapidities. The change in behavior is due to large dipole configurations, 
which has been discussed above. To illustrate the behavior of A^ as a function of the impact 
parameter, we have parameterized our solution as: 

where S{b) = {R^ + and / and g are fitted functions. Eq. (fTUI) is presented in Fig.Elas a 
dashed line. In our approach, since we exclude large dipoles from the evolution, the amplitude 
at large y deceases with b faster than 6""^. As was pointed out in Ref.jH] the power-like decrease 
leads to power-like growth of the average 6^ as a function of energy. However, Fig. [7| shows that 
at large values of the dipole size there is only a slow increase, (if at all). Indeed, the simple 
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Figure 4: The evolution of as a function of r/R, for fixed b. 

formula (6^) = a + hy + cy"^ provides a good description of the enegy behaviour of the radius 
of interaction, for different values of r. (see Fig.[7j). 

We have also used Eq.()10|) to parameterize the correction AA^^^-* = / ANj^^'^dy. In Fig.|Hlwe 
show AA^'-^-* and its parameterization for large r/R = A and large y = 14, we stress that 
the power tail originates from the correction to the first iteration. 

Fig. El and Fig.^J demontstrate the essential difference between running and fixed QCD 
coupling cases. One sees that for running QCD coupling the amplitude has a steeper growth 
with energy, reaching the saturation boundary for smaller y than in the case of fixed as- On 
the other hand, the behaviour of the tail in the impact parameter distribution looks very 
similar in both cases (see Fig.[TT|). It should be stressed that the scattering amplitude falls off 
exponentially at large b { N oc e~^^ with m = 1 to 2 GeV, see Fig.lTTl Fig.|Sland Fig.lTTl 
show that our solution has no power-like tails in the h distribution, and therefore, supports 
our hypothesis that the power-like large h decrease, stems from the dipole configuration which 
cannot be treated within the non-linear Balitsky - Kovchegov equation. 



3.2 Saturation scale 

We shall now use our solution to the non-linear evolution equation to investigate saturation 
effects and scaling properties of the dipole- nucleon scattering amplitude. First, we calculated 
the saturation scale r^d/, h) as being the value of the dipole size, at which the amplitude 

Ar(r„ 6) = 1 - e-^ = 0.632 (11) 

Fig. IT^ shows the solution to Eq. (fTT|) . We recall, that for the case of fixed QCD coupling in 
the semi classical solution (see Ref. pilj) we expect that at large y and b 

Ql{y,b) oc e'ye-^, (12) 
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b = 




Figure 5: The energy dependence of at 6 = 0, for different values of r/R. 

where the estimate for the value of A, is A = Aas and R is the radius of the target. 

For the running QCD couphng constant, the situation is more interesting (see Ref. [3 EHl 
113]): (i) as function of y 

q2 ^ (13) 
in the region of b < RXri/ ; and (ii) for 6 > RXru 

Qs oc e"^ ; (14) 

with \r ^ 10. 

To aid in understanding the main features of the energy and impact parameter dependence 
of the saturation scale, we fitted the result for the saturation momentum using the following 
parameterizations, which reflect the above semiclassical properties: 

At fixed 6: Qs{y,b) = Cy^ie^^^ + C^,2 eV^^^ (15) 
At fixed y: Q'Ay,b)/4: = ^^sh^ (16) 

In Fig. El we show how these parameterizations describe the calculated saturation momen- 
tum for fixed b. For running QCD coupling we used the value of A/j = 0.26 to 0.30, which is 



9 




Figure 6: The impact parameter dependence of for r/R=l. The dashed line corresponds 
to the parameterization of Eq. ((TU)) . 

much smaller than the semiclassical prediction: Xr = 10. For fixed QCD coupling, we see an 
exponential growth with the value of Ai? ~ 1.4 to 1.6, which is approximately two times larger 
than the theoretical prediction, Xp = Aas ~ 0.8. 

Both Xp and Xr turn out to be independent of b. 

Fig. IT^ shows the impact parameter behavior of the saturation scale. Both y and b behaviour 
of the saturation scale are in a qualitative agreement with the expectations mentioned above. 
The formula given in Eq. ()16|) describes the main features of the 6-behavior rather well. For 
fixed QCD coupling, we see that for a wide range of y we have a Gaussian behavior, while 
for the case of running QCD coupling the decrease sets in, only in the region of very large b. 
However, the surprising fact is that even at fixed as for large y, the 6-behavior tends to the 
same distribution as in the running as case. This observation requires a qualitative explanation. 
Note that at large b in both cases, the impact parameter dependence appears to be Gaussian, 
with a radius of interaction which increases with energy. In our parametrization the radius is 
equal to RinteracUon = C3 R and C3 (X y. 
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Figure 7: Energy behavior of the average h at different values of the dipole size for fixed 
(sohd) and running (dashed) QCD coupling. 



In Fig.^Jthe large h behaviour of the saturation scale is plotted. Qs falls off as e~°'^'' . One 
expects such an exponential behaviour from general properties of analyticity and crossing sym- 
metry 1 and this behaviour supports our point of view that the initial conditions determine 
the large h behaviour of the scattering amplitude. 



3.3 Geometrical scaling 

In |21 m 1^ it was shown that the dipole-nucleon scattering amplitude exhibits scaling, which 
occurs due to the relationship between the y and r dependences of iV. This scaling behavior 
was confirmed in the approximate solution of the non-linear equation [21], without theoretical 
knowledge of the impact parameter behavior. More specifically, the scaling property implies 
that, in the kinematical region where saturation occurs, the amplitude is a function of one 
variable r = rQsiy, b), with Qs = 2r7^ 

In Fig.Hniour solution is plotted as function of r. A glance at this figure shows that the 
solution displays the scaling properties at large values of y, as was predicted by theory j2I]. On 
the other hand. Fig. clearly displays that geometrical scaling behavior is more pronounced, 
for the running QCD coupling constant case. 

The improved geometrical scaling which we observe for the running coupling may, at first 
sight, be confusing, since contains a new dimensional parameter Aqcd, which in principle 
could spoil the scaling. However, in Ref. |13j it has been suggested that the running of the 
QCD coupling freezes at the saturation scale, thus recovering the broken geometrical scaling 
property. We interpret our numerical solution, as supporting this hypothesis. 

Nevertheless, as the exact definition of the saturation region is a debatable point, we shall 
also investigate the scaling property of N, using the ratio: 

^_ dN{y,r,b)/dy 
d\n N{y, r, b)/dr 

If scaling does exist, TZ should equal d\nQs{y,b)/dy. Note, that the scaling property of 
implies that in the saturation region, TZ should not depend on r. 
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Figure 8: The impact parameter dependence of AA^. Tlie dashed line corresponds to the 
b~'^ parameterization of Eq. (|1U() . employed for AA^^^^. 

In Fig.Elwe present the ratio 7^ as a function r, for different values of y and b. Although, 
as stated, our solution is only applicable at small values of r, we see that for long distances (for 
r/R> 0.07 with R = \/3.1 GeV~^)) TZ is reasonably flat as a function of the dipole size, while 
for smaller values of r the ratio falls steeply as a function of r. Once more, the geometrical 
scaling behaviour is more striking for the running QCD coupling case. The second interesting 
observation is that this behaviour is concentrated at small b, while at large values of b we see 
a strong violation of geometrical scaling. 

3.4 Comparison with other solutions 

The solution in j6] was the first step at solving the equation for fixed impact parameter. This 
solution provided an excellent reproduction of the F2 HERA data with x^/n.d.f. = 1, using 
i?2 = 3.1 GeVl In Fig.dl we compare our solution (with the initial conditions of Fig.|2b) and 
the solution of ^ (with the initial conditions of Fig.|2t)- Fig. 118b. shows the amplitude as a 
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function of r / R. The r-dependence of the solutions is substantially different. We see that the 
6-dependence of the equation generates a rapid rise of the amplitude as a function of r. 

Fig.lTHb shows the 6-dependence of the two solutions. The impact parameter dependence of 
our present solution is the result of the evolution process (starting from the input of Fig.|2l3), 
whereas the impact parameter dependence of was put by hand, according to the Glauber- 
Mueller formula N{h) = 1 — exp(— kS'(6)), with n = — ln(l — N{b = 0)). The generated 
6-dependence differs from the ansatz of 6J. For low energy {e.g., y = 5.6) the fixed b solution 
extends to larger values of the impact parameter, whereas at large energies {y = 9.8) our 
present solution lies well above the fixed b solution for values of b about up to 6 = 10, where 
both solutions are small. Only for intermediate energies and distances {y ^ 7, r ^ -R/2), do 
the two approaches exhibit similar behavior. 

In Fig.lTHh we show how the two solutions differ in terms of the total cross section, i.e., the 
area under the curves of Fig.lTHb. Only for y ~ 7 and small distances do we find similarities 
between our present solution and the fixed b solution of [Hj. 

Fig. IT^ shows a qualitative comparison of our solution and the solution of Ref. ^1]. Only a 
qualitative comparison can be made here, as (i) in |14j there was an explicit cos 6 dependence, 
where 9 is the angle between r and b, whereas our solution is an integrated quantity, N{y, r; b) = 
/ df\\N{y, r; 6; f\\); and (ii) there is a normalization difference, due to the different scales at which 
the dimensional variables are defined. 

The origin of the last difference is the different initial conditions that were used. In 
a pure Glauber-Mueller form for the initial conditions was assumed, using a profile S{b) = 
lOexp (—6^/2). Our initial conditions, on the other hand, have a different shape (see Fig.|2l3) 
and are based on parameters which were established for our fits to the HERA data jH]. 

To enable us to make an illustrative comparison, we have used the following procedure: 

1. To minimize the angular sensitivity, we compare the two solutions only at low r or low b. 

2. We interpreted the input impact parameter dependence of [Hj ( i.e. S{b) oc exp(— 6^/2) ), 
as exp —b'^/R^ with R^ = 2. Recalling that we use R^ = 3.1 GeV^ as the size of the target, 
we have rescaled b according to 6^ — > 26/3, hence matching the dimensional units of 

3. The shape of our initial input (see Fig.|2l3) is substantially different from the Gaussian 
of On the other hand, we wish to compare the solution at y > yo, where we assume 
that has a weak dependence (if any) on the shape of the initial conditions. Thus, we 
have chosen to rescale r so that after evolution of, say, two units of rapidity^, the two 
solution will match at small b and fixed r. We found that a rescaling of r — > 2r is 
sufficient to produce such matching. 

^Note that there is also a difference in the definition of y: in yo = and here yo = — In 10^^ « 4.46. 
The values of y in Fig. ^1 correspond to the latter definition 
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Referring again to Fig. we see that as a function of b, our solution decreases more rapidly. 
In addition, we see that as a function of r and 6 = 0, the large normalization difference between 
the two solutions at y—yo = 2 diminishes at sufficiently large y—yo = 9. The 6-dependence of the 
saturation scale Qg is similar, but there is a considerably large normalization difference, which 
may be due to disparities in the definition of the saturation region, and the above mentioned 
dimensional parameters. 

4 Summary 

In this section we summarize the main features of our numerical solution of the non-linear 
evolution equation. 

• It is shown that the non-linear equation in the region where it's validity is guaranteed, 
leads to a solution which falls off as 



• We showed that the inclusion of dipoles of size larger than the target (R), generate 
the power-like tail in impact parameter distributions. This means that the power-like 
behaviour is an artifact of using the non-linear equation, for calculating the scattering 
amplitude in the kinematic region, where this equation is not valid ( see Refs. [23] for an 
approach that can describe this region); 

• In our approach we modified the kernel of the non-linear equation by introducing Q{R— 
Ti) Q{R — r2). We concur with the idea of Ref.j^l that to obtain the correct (exponentially 
decreasing) large impact parameter dependence, it is necessary to change the kernel of the 
non-linear equation. On the other hand, the non-linear equation was obtained implicitly 
assuming that the size of the projectile dipole is much smaller than the size of the target. 
This was a misconception based on the semiclassical approach, and on the exact solution 
for the simplified double logarithmic kernel for r -C i?, that large dipoles will not appear 
in the intermediate stage of evolution ^TJ |Q. We demonstrated that the large b tail 
which behaves as 1/6^, stems entirely from the large dipole contribution, even for small r. 
The large b tail gives a very small contribution to the amplitude. One can compare Fig.lTUl 
with the value of AA^ in Fig.|Hl to realize that AN/N is negligible for 6 > 25 GeV~^. 
Therefore, the statement of Refs. P^IT^ looks reasonable, that the large b tail which is a 
problem in the pQCD approach, does not contributute to the main physical observables, 
including the total cross section. 

• The value, energy and impact parameter dependence of the saturation momentum are 
calculated. We find that at large values of rapidity, the value of the saturation momentum 
no longer depends on the impact parameter b. A behaviour of this type was predicted for 
Qs in the case of running QCD coupling [2H], but was not expected for fixed as', 



N{r, b;y) oc e r for large b > Xy ; 
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• The geometrical scaling behaviour of the scattering amplitude is discussed, and we show 
that the scattering amplitude can be displayed as a function of one variable r = r Qsiy, h). 
The geometrical scaling behaviour is more striking for the running QCD coupling, than 
for fixed a^. An explanation for such behaviour could be related to the idea ^21; that 
the running QCD coupling is frozen at the saturation scale. In this case the smallness of 
asiQs) makes our calculations more reliable theoretically, and leads to better geometrical 
scaling; 

• The radius of interaction increases with energy, approximately behaving like y < b > 
(X y, in accordance with the Froissart theorem [lOj. This result confirms that the large b 
behaviour stems from the non-perturbative scale introduced in the initial conditions; 

• The solution that we found is quite different from previous attempts, where an anzatz 
for b dependence was assumed. We plan to use our solution for developing a new procedure 
for a global fit to DIS data; 

• We stress that the solutions for running and fixed as are quite different. For running as 
we obtain a larger interaction radius (approximately twice as large), steeper y dependence, 
and a larger value of the saturation scale; 
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Figure 9: Our solution for the scattering amplitude as function of the impactparameter b, 
at different energies and at different values of r for running QCD coupling. 
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Figure 10: Our solution for the scattering amplitude N as function of the impact parameter b, 
at different energies and at different values of r for fixed QCD coupling. 
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Figure 12: The saturation scale as a solution to Eq. ()11|) of the amplitude for fixed and running 
QCD coupling, for different values of y and h. 
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Fixed tts 



Running a. 




Figure 13: The energy behavior of the saturation momentum Qg — 2/rs for fixed and running 
QCD coupUng. For small y the curve is InQs = a + ^ |/, while for large values of y the curve 
corresponds to In = a + | ^/\r+T). The results of our calculations are presented as circles 
whose magnitude reflect the size of possible numerical errors. 
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Figure 14: The impact parameter (6) behaviour of the saturation momentum Qs = 2/rs for fixed 
and running QCD couphng. The curve correspond to Eq. (fT5|) .The results of our calculations 
are presented as circles whose magnitude reflect the size of possible numerical errors. 
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Figure 15: Large impact parameter (6) behaviour of the saturation momentum Qg for fixed and 
running QCD couphng. The curve corresponds toQg — C e~'^^ . The results of our calculations 
are presented as circles whose magnitude reflect the size of possible numerical errors. 
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Points: R= (dN/dy) / (dN/d(lr r)] after smoothing, b=0 . Points: R= (dN/dy) / IdN/d (In r] ) after smoothing, b=0. 
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Figure 17: Scaling ratio 7?. as a function of dipole size r for diflFerent values of y and h. 
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Figure 18: Comparison of our solution (solid) and the fixed b solution [2j (dashed). The com- 
parison is for (a) the r-dependence of the amplitude at 6 = 0; (b) the 6-dependence 
of the amplitude at r = R/2; and (c) the integrated amplitude, / (PbN{y,r,b) as 
a function r/R. 
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Figure 19: Comparison between our solution (solid) and the solution of 14^ (dashed). 
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